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ABSTRACT 


This thesis presents a new approaeh to solving a elass of problems arising in the 
design of satellite swarms. Using the fundamentals of optimal eontrol theory, a 
framework is developed that eaptures the essenee of “eoneurrent” design and eontrol of 
spaeeeraft formations. This framework is used to artieulate a variety of formations 
ineluding the notion of an aperiodie formation. Additionally, formations that require 
aetive eontrol are presented along with their eorresponding thrust profile. Based on a 
deliberate problem formulation, whieh ineludes mass as a state variable, it is shown that 
the numerieal approaeh easily handles nonlinearities. Using the general-purpose dynamie 
optimization software, DIDO, this thesis demonstrates how a minimum-propellant 
formation eonfiguration ean be easily designed for satellite swarms without the use of any 
analytical results. If a zero-propellant eonfiguration does not exist, then this method 
automatieally determines the minimum fuel and the assoeiated eontrols required to 
maintain the eonfiguration. This thesis lends eredenee to the notion of numerieally 
searehing for minimum-fuel formation eonfigurations for spaeeeraft swarms subjeet to 
arbitrary nonlinear dynamies. Thus, praetieal formations may be designed and 
eontrolled using this method. 
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I. INTRODUCTION 


The concept of distributing the functionality of larger satellites among smaller, 
cooperative satellites has been seriously considered for assorted space missions to 
accomplish goals that are not possible or very difficult to do with a single satellite. ’ ’ 
The current trend in satellite design is toward using smaller, more capable satellites in 
cooperative formations or distributed arrays."^’^ 

An extremely critical issue in architecting satellite swarms is the formation 
configuration. Arguably, the main feasibility criterion in the architecture of satellite 
swarms is the design of globally-minimum-fuel configurations. This simply follows from 
the fact that fuel for orbiting spacecraft comes at a very high premium and could 
significantly offset any other advantage held by a swarm configuration. Thus, there is a 
need to determine zero-propellant formation configurations (if they exist) and methods 
for controlling the formation with little or no propellant. It is well known that a family of 
zero-propellant circular and elliptic formations exists when the spacecraft are subject 
only to an inverse-square gravity field. ’ However, these formations tear apart in the 
presence of “disturbing” effects such as J 2 . Thus, a search for invariant relative orbits or 
formations (if they exist) goes on. Another effect that must be accounted for is a non¬ 
zero eccentricity of the reference orbit.^ Unlike the J 2 disturbance, an error in eccentricity 
can be controlled by a one-time expenditure of propellant. This is based on the 
observation that the J 2 disturbance is an error in the dynamical model whereas the error in 
eccentricity is one of initial conditions. The error in eccentricity does not fully address 
the problem since in many applications it is desirable to have formations for every 
eccentricity (and not just small eccentricities). Hence, the “real” problem is to find 
formations in the presence of the totality of (modelable) deterministic forces. If such 
formations do not exist naturally, then it is imperative that the minimum-fuel formation 
configuration be determined. The total fuel consumption of a formation over its lifetime 
then establishes the practical feasibility of a swarm architecture. There are some 
opinions^® suggesting that minimum-fuel configurations exist and that they do not offset 
the advantages offered by other performance metrics. However, there appears to be no 
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general-purpose method published in the open literature on how to find these minimum- 
fuel eonfigurations subject to arbitrary forces. 

This thesis presents a general problem algorithm for finding optimal formations 
and shows that a very natural setting for solving this problem is (nonlinear) optimal 
control theory. In particular, these ideas are formulated by using elements from optimal 
periodic control theory with partially periodic states}^ The Legendre pseudospectral 
technique is applied to numerically solve this problem (see Reference 12 and the 
references contained therein). Because of the Covector Mapping Theorem, this 
technique is neither a direct nor an indirect method. Rather, it provides all the ease of a 
direct method while providing the accuracy of an indirect method. This method is 
implemented using the general-purpose software package, DIDO,'"^ which has been used 
extensively over the past few years to solve a myriad of complex optimal control 
problems. 

The applicability of optimal control theory to satellite formation was not known at 
the outset of the thesis work presented here. For this reason, the research style chosen 
was a building-block approach. First, optimal control theory had to be demonstrated as 
an appropriate framework for a simplified model and then complexity would be added 
incrementally. 

A model simply captures the essential aspects from a certain point of view and 
simplifies or ignores the rest. The so-called Hill-Clohessy-Wiltshire (C-W) equations 
were chosen as the first model specifically because the solutions were known. This 
allowed a validation of the method before embarking upon more complicated models. 
The second model was chosen to address the eccentric reference orbit. To do this, a 
nonlinear version of the linearized equationswas used for the dynamical model 
without the presumption of a solution. 
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II. GENERAL FRAMEWORK 


A. PROBLEM DEFINITIONS^ 

The notion of a swarm can be defined by stating that a group of satellites are said 
to be in formation if a given configuration metric c(|x'|) is bounded, 

c, <c({x'})<c„ (2.1) 

where x' denotes the state of the spacecraft at time x. This state can be the usual 
position-velocity state or any other set (e.g. orbital elements). 

The dynamics of a swarm of Ns spacecraft are given in some coordinate system 

by, 

x'=f(x>',x;p) i = \...N^ (2.2) 

If the distance between any two satellites, J(x',x^), is chosen as the metric, then 
the swarm is said to be in formation if 

df <d{x‘,x^')<d‘f Vx (2.3) 

where di and du define the smallest and largest allowable separation distances, 
respectively. Instead of choosing separation distances between every spacecraft pair, it is 
sometimes simpler to choose a separation distance between a spacecraft and a reference 
spacecraft. In this case, equation (2.3) can be replaced by 

dj <d{y,x^)<di Vx (2.4) 

where y is the state of the reference spacecraft. From these fundamentals, a family of 
formations can be defined as follows. If d{y,x^) is a constant for all j, then the 
formation is a circular formation, 

d{y,x^) = dl =di Vx (2.5) 

with the spacecraft at reference point y called the “mother” and the remaining j 

spacecraft denoted as “daughters”. A circular formation can be defined even in the 

absence of a mother spacecraft. Thus, the mother spacecraft may be replaced by a 
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reference point, y, which serves the purpose of providing a non-inertial reference frame 
to the entire formation. Multiple rings of circular formations can similarly be defined 
with multiple distance bounds on a collection of spacecraft with respect to a single 
reference point. 

A formation is defined to be fully periodic if the entire state vector is periodic, 

x'(Xo) = x'(x,) (2.6) 

and partially periodic if only some of the components of the state vector are periodic. 
For example, a formation may be partially periodic because it is periodic in position but 
not in velocity, or vice versa. If propellant is used to maintain a formation, then by this 
definition, the formation is partially periodic if the aperiodic mass is included as a state 
variable. It is apparent that a circular formation is a periodic formation but the reverse is 
not necessarily true. Of special note is that this definition for periodic motion can be 
either inertial or relative. Hence, periodicity in the relative frame does not necessarily 
imply periodicity in the inertial frame. That is, the swarm may drift in inertial space, but 
it will stay cohesive as a formation. 

For a fully periodic formation using Cartesian position and velocity vectors as the 
state, the periodicity constraint may be written as, 

r(Xo) = r(x^) (2.7) 

f(Xo) = f(x,) (2.8) 

These conditions allows us to further define two classes of partially periodic problems: 

(1) when only equation (2.7) is imposed while the boundary conditions on f are free and 

(2) when only equation (2.8) is imposed while the boundary conditions on r are free. 
This thesis will limit its attention to fully periodic solutions and the two classes of 
partially periodic formations described above. 
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Finally, it must be noted that a formation need not be periodie at all! Henee, a 
relaxed formation is defined as the ease when the state veetor returns to within a defined 
spaee around the initial state after some time , 


f < X* — x" < f 

- Aq a^ _ 


Thus, the familiar notion of an epsilon ball. 


x' -x' 

0 f 


<£ 


(2.9) 


( 2 . 10 ) 


is ineluded in this definition. 


The eonfiguration is eonsidered to be optimal if, in addition to satisfying the 
eonfiguration constraint, a scalar performance measure. 


1 ' 

J[x(-), u(-), To, x^; p] = —— J F(x(x), u(x), x; p) Jx 


( 2 . 11 ) 


is minimal. The reason for choosing a cost functional borrowed from optimal periodic 
control theory is that orbital motion is, by nature, periodic. Further, in addition to finding 

minimal fuel configurations, it is also desirable to find the optimal period, [x^ that 

renders a minimal cost per period. 

To complete the optimal periodic control formalism,the configuration 
constraints described in equation (2.1) are broken down into event constraints and path 
constraints. Event constraints are end point boundaries defined by 

e, ^e(x'(Xo),x'(x^),Xo,x^)<e„ (2.12) 


Path constraints are boundaries placed on the trajectory of the model, 

g/^g(x' (x),u' (x),x)<g„ (2.13) 


Additionally, each of the state and control variables may have a constraint placed on 
them by 
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( 2 . 14 ) 


X* < x” (x) < X* 

u‘<u'(x)<u' (2.15) 

All of the constraints shown above can be used as equality constraints by simply 
setting the upper and lower bounds equal. They are written as inequalities for the 
purpose of generality. Any formation configuration may now be defined in this 

“standard” form as finding the controls, u'(x), and the optimal period, [xy^-Xg], that 

minimize the cost of equation (2.11). 

1. Reference Frame 

In order to describe relative position, motion, and configurations, the Formation 
Reference Frame will be used throughout this thesis. Figure II-l shows this reference 
frame, which is defined with x pointing in the radial direction, y pointing perpendicular 

to X along the direction of motion, and z completing the right-handed coordinate 

20 

system. This reference frame is often called RSW or Satellite Coordinate System. 
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The assumption of the eoordinate system is important in defining any problem. 
For some eoordinate systems, the formation relative motion and eonliguration eonstraints 
are intuitive. For example, the path eonstraint for a eireular formation ean be deseribed 
in the formation referenee frame simply using the relative position 

+ r'^ = constant (2.16) 

On the other hand, to deseribe this eonstraint in the inertial frame requires a 
transformation to that frame. This transformation need only be eompleted onee, but ean 
be very labor intensive. 

This same eoneem must be applied to the dynamie equations of motion. The 
eomplete, non-linear equations ean be readily expressed, with no assumptions, in the 
inertial frame by 

(2.17) 

R. 

Conversely, in the formation referenee frame a transformation is again required and 
usually involves introdueing assumptions and linearizations. 

Equal to the eoordinate system in importanee is the ehoiee of the variables used to 
deseribe the satellite state. One set of variables is the Cartesian position and veloeity 
veetors. Other sets available inelude many different orbital element sets, whieh are 
espeeially useful if the eoordinate system is Earth eentered and inertial. However, if the 
solution ealls for eontrol thrusting, it will be a eomplex translation into the orbital 
elements. Eor these and other reasons, all models presented in this thesis utilize Cartesian 
position and veloeity veetors in the Eormation Referenee Erame as the basie spaeeeraft 
state. Depending on the model, there may be additional variables in the state, but there 
will be a minimum of these six. 

B, SOLVING OPTIMAL CONTROL PROBLEMS 

Until reeently, solving general nonlinear optimal eontrol problems was an arduous 
or impossible task. The theoretieal framework for solving sueh problems is the Minimum 
Prineiple. Numerieal methods based on the Minimum Prineiple are known as indireet 
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methods. While solutions obtained from indireet methods are aeeurate in the sense that 
they satisfy the neeessary eonditions of optimality, they are fundamentally burdened by 
numerieal sensitivities as noted by Kalman over four deeades ago. The so-ealled 
indireet multiple shooting methods and indireet eolloeation methods overeome this 
eomputational instability problem but at the expense of eonvergenee: good guesses on the 
eostate time-history are neeessary to suoeessfully solve the problem. Over the last 
deeade, the so-ealled direet methods have eome to the forefront. These methods bypass 
the Minimum Prineiple and “direetly” solve the problem in various ways. Betts provides 
an exeellent review of this approaeh. Early direet methods were plagued by 
inaeeuraeies, partieularly in the determination of the eontrols. More reeently, major 
breakthroughs in higher-order methods and large sparse numerieal methods have quiekly 
narrowed this gap. One partieular approaeh is to use a solution obtained from a direet 
method as a guess for an indireet method. Another approaeh, favored in this thesis, is 
ealled the Legendre pseudospeetral method. This method is used to solve the formation 
design and eontrol problems posed and is implemented in the reusable software paekage, 
DIDO. Unless otherwise speeified, all results reported in this thesis are obtained using 
this software. 

1. Necessary Conditions for Optimality* 

The first step in solving an optimal eontrol problem is to eonstruet a sealar 
funetion ealled the Hamiltonian, H , 

//(x,u,t,k) = —f(x,u,t) ( 2 . 18 ) 

X f — Tq 

where f(x,u,t) are the dynamie eonstraints on the system, and k(t) are the Lagrange 

multipliers ealled eostates. Aeeording to the Minimum Prineiple, at eaeh instant of time, 
the optimal eontrol is obtained by solving the following problem. 

Minimize H with respeet to u , with u e U 


* Most of the information in this section comes from class notes and discussion from AA 4850 and is 
reproduced here for completeness. 
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where U is the set of all allowable eontrols. To solve this problem, the Lagrangian of 
the Hamiltonian must be eonstructed; 

T(x,u,t,J^,(|)) = //(x,u,t,^) + (t)(t)^g(x,u,t) (2.19) 

where g(x,u,t) are the path eonstraints and are the assoeiated Lagrange 

multipliers. The path eonstraints inelude all trajectory path constraints as well as any 
state and control bounds. Applying the Karush-Kuhn-Tucker (KKT) theorem to the 
Lagrangian results in a set of necessary conditions and provides a means to demonstrate 
the optimality of a solution. 




du 

(2.20) 



(|)(t)^g(x,u,t) = 0 

(2.21) 

with 

<0 

g, =g(x,u,x) 


(]) = < 

if 

= 0 

g(x,u,x) = g„ 

term by term 

g, <g(x,u,x)<g„ 

(2.22) 


any 

gl=gu 


The third case above describes a 

special condition when the constraints in g 

are interior 

or non-binding. 


g/<g(x,u,t)<g„ 

(2.23) 

For these cases, the multipliers, (|) 

= 0 and equation (2.20) simplifies to 




dL _ a// _ ^ 

du du 

(2.24) 


It is desirable to have interior constraints because the problem behaves as if the 
constraints do not exist. For this reason, the constraints placed on a problem may have an 
actual value in practice, but if sufficiently large as to be non-binding, they can be 
described as unconstrained. 
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2. Scaling 

Scaling is critical to all optimal control problems. The goal is to establish a 
seheme that seales all parameters and state variables so that they are elose to one another 
in value. For example, it is better to have positions from 0-10 with veloeities from 1-10 
than positions from 7000-17000 with veloeities from 0-10. The method for sealing the 
problem starts with setting a desired unit and determining the remaining units that 
eomply with this standard. For the problems that follow, the desired standard unit was 
time. The time unit was ehosen to go from 0 to 1 for a single orbital period. The 
normalization eonstant, TU , is defined to be equal to the period of the referenee orbit in 
seeonds. This ean be done several ways, but the most simple way is to seleet a value for 
the semi-major axis, a , and use the following equation 


TU = 2n 



(2.25) 


where |l is the gravitational parameter of the earth and is set to 


11 = 3.986004418x10'" ^ 

s 

The normalized time now beeomes 


(2.26) 


T 


nondim 


T 

TU 


(2.27) 


and therefore, the orbital period beeomes 1.0 TU. This definition for the time unit, by 
nature, defines the mean motion sinee 


n 


nondim 


2k 

Orbital Period 



(2.28) 


The next variable to be seleeted is the mass unit, MU , whieh is simply ehosen to 
be equal to the initial mass of the spaeeeraft or 




-^ = 1.0 


MU 


(2.29) 
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The last variable that is ehosen by the user is the standard distanee unit, RU . The 
ehoiee for this unit is based on the desired formation spaeing. It ean be varied to provide 
larger or smaller formations without affeeting the eonfiguration. Given TU and RU , 
the remaining variables ean be normalized as follows 


^ nondim 


r 

RU 


(2.30) 


nondim 



(2.31) 


a 


Pmmdim 



(2.32) 


Unless dimensional units are speeified, all values are assumed to be sealed and 
nondimensional and the nondim subseript will be dropped. 

3, Unique Issues for Periodic Problems 

As with any optimization problem, most errors are not in the solution but in 
asking the right question. One speeial issue with periodie optimization problems is 
formulating an initial guess. Every optimizer requires a first guess, whether it is provided 
by the user or ealeulated by the solver. The quality of the guess is often an issue, as some 
solvers require relatively aeeurate guess. DIDO does not require a good quality guess or 
even a feasible guess. However, a simple linear interpolation between the endpoints is 
not an option for periodie problems. The explanation lies in the periodieity eondition 
itself. For orbits, even a straight line eannot be used as a guess sinee the motion is 
elliptioal in nature. Some sort of ellipse must be used as a guess for orbital motion. 

Determining the eost funetion is another diffieult part of eonfiguration design. It 
affeets not only the speed of the solution, but the solution itself For example, beginning 
with seetion III.D, thrust and mass are used in the non-linear equations of motion. The 
eost funetion ean be written simply as 


J 


Mj- —Mq 


(2.33) 
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On the other hand, it ean also be written as 




Xf—X 

f 0 -to 


(2.34) 


whieh also results in a minimum mass expended. It turns out that for zero or very low 
thrust problems the ehange in mass is very small and may be below the numerie toleranee 
of the solver. This often results in a sub-optimal solution. When thrust is used in the 
eost, sinee it is a eontrol variable, the solver is better able to find the optimal solution. 

Another aspeet that proved interesting was the -Tq) term, used to represent 

the optimal period. For Natural or Zero-Thrust solutions, the optimal period was exaetly 
equal to the orbital period. However, onee thrust is required to maintain the formation 
the optimal period may or may not equal the orbital period. One way to get around this 
disparity is to foree the solution period to be equal to the orbital period, that is =1.0. 


C. VALIDATING SOLUTIONS 

The purpose of this seetion is to deseribe the various methods used to validate the 
numerieal solutions. Numerieal methods may seem to diseard the physies of the problem 
and rely solely on the mathematies, but this is far from the truth. Any solution found 
numerieally must be ‘assessed’ to see if it is in faet a feasible and legitimate solution. It 
is preeisely here that the “Physies” intuition and knowledge are implemented to assist in 
verifieation. To that end, there are many different ways to verify that a solution is 
feasible. The optimality of a solution ean be demonstrated using the optimality 
eonditions deseribed in seetion B.l. Eaeh of the three primary means of validating 
solutions are deseribed below. 

1. Numeric Propagators 

The primary method for validating feasibility of the solutions is to propagate them 
using a numerieal propagator. Sinee the equations of motion are Ordinary Differential 
Equations (ODE), a tool is needed to solve them for eaeh desired time step. The 
environment for DIDO is MATLAB®, so a resident ODE solver was used. There are 
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several to ehoose from, but ODE45 is the most oommon and usually the most eapable for 
non-stiff equations. 

The first step in validating is to build a propagator from the ODE solver. This 
requires defining the Equations of Motion (EOM) that govern the behavior for the solver 
and establishing the initial eonditions. The initial eonditions were taken direetly from the 
DIDO solution. It is important for the EOM to be the same as the ones used to find the 
solution, to ensure an aeeurate validation under the same assumptions. It does no good to 
make assumptions in defining the problem, only to forget about them in validating the 
results. 

2, Commercial Software 

Another method employed to validate solutions is the use of oommereial software. 
Both Satellite Tool Kit® (STK), from Analytieal Graphies, Ine. and EreeElyer®, from A.I. 
Solutions were used as an independent souree for propagating the formations. Usually 
these programs are used for visualization and presentation, but both software suites 
inelude a very robust non-linear propagator. Any number of reports ean be seleeted to 
identify and traek the desired parameters. Eaeh of the programs also ineludes a viewing 
eapability that allows the user to be plaeed on the referenee point observing the formation 
from that perspeetive. This direetly eorresponds to the Eormation Referenee Erame 
shown in Eigure II-1 and is very useful in visualizing relative formations. 

The meehanism for validating solutions using these programs is to import the 
initial conditions for the individual formation satellites into the programs and propagate 
them forward. This ean be done direetly from MATLAB for use in EreeElyer, or with an 
Ephemeris file for use in STK. If the solution requires aetive eontrol, then the eontrol 
history must also be imported into the programs. Other than initial eonditions and/or 
eontrols, there is no other information provided to the software. This laek of information 
is exaetly what validates the solutions. If the formation behaves as predieted in the 
solution when propagated by the software, then the formation is feasible. 
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3. Previously Discovered Solutions 

The third method is used to validate the method more than it is to validate a 
partieular solution. If the method is employed and is able to identify previously 
diseovered solutions, it is valid. The ultimate validation would be to demonstrate that 
this method is able to find all previously diseovered solutions. Rather than devote 
preeious time to an exhaustive eatalog of all prior solutions to all satellite formations, 
several representative formations were ehosen as demonstrations of the ability. 
Speeifically, the well-known eireular and projeeted eireular solutions to the Hill- 
Clohessy-Wiltshire equations were reprodueed. For an example using elliptieal referenee 
orbits, Referenee 17 outlines a solution with e = 0.7 for the referenee satellite. 
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III. CIRCULAR REFERENCE ORBIT 


A. PROBLEM DEFINITION 

The framework deseribed in the previous ehapter ean be applied to designing and 
eontrolling spaeeeraft formations subjeet to any arbitrary forees. This chapter begins by 
applying this framework to spacecraft subject to an inverse-square gravity field only. If 
the equations of motion are linearized in the formation reference frame about a circular 
reference orbit, the C-W equations are obtained and closed-form solutions are easily 
found. From these equations, it is apparent that zero-propellant formations exist making 
this problem an excellent starting point. These equations also served as a good model 
with which to validate the process for this class of periodic problems. 

One zero-propellant solution to the C-W equations is a circular formation. In 
order to solve this problem the following configuration was used. The state consists of 
the Cartesian relative position and velocity components shown below 

x = [r' ^ K K ^7 (3.1) 

The state constraints were 


-2r < r r r < 2r 


-r <. r r r < r 
max — a: ’ v ’ 2 — max 


(3.2) 


where was chosen arbitrarily. The position constraints, while specified, were never 
active due to a tighter path constraint. The velocity constraints were not active either due 
to the choice of . The controls, representing the relative accelerations in each axis, 
were constrained by 


—u <u ,u ,u <u 

max — 7 — o 


(3.3) 


with specified by the user. Since the solution is a zero control solution, these 
constraints were also inactive. 
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The event eonstraints were specified for a fully periodic solution using equations 
(2.7) and (2.8) although they were written in the form 

x(To)-x(t^) = 0 (3.4) 

or more specifically 

fi(v) 

t(y) 

The path constraint, g , for a circular formation is defined by 

{r-bf <[r^+ry+r^^)<{r + 8f (3.6) 

This general form of the path constraint includes an allowance for a tolerance, which for 
circular formations is set to zero or for near-circular formations is nonzero. For this case, 
the path constraint for a circular formation and the event constraint for a fully periodic 
formation are redundant. The path constraint will force a fully periodic solution, even in 
the absence of any event constraints. This was demonstrated and validated during several 
of the solution sets. 

The cost function to be minimized was 

J =- {(u^ + + u^ \dx (3-V) 

X , -X ^ ' 

f 0 -to 

which is identical to a cost that was based on the absolute value of each control 
acceleration, if the solution is a zero-control solution. 



"fiK)" 

My). 
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B, EQUATIONS OF MOTION 

Since the equations of motion for this model are well known'their 
derivation will not be described here. Instead, only the equations and their assumptions 
will be presented. 



■ 0 

2n 

0" 


~3n^ 
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0 " 



f = 

-2n 

0 

0 

f-1- 
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r-i- 

u, 


0 

0 

OJ 


0 

0 

2 

—n 


_“z_ 


(3.8) 


where r is the vector representing the inter-satellite distanee expressed in the formation 
reference frame (see Figure II-l), 


r = [x y zf 

or in state form: 


' 0 

0 

0 

1 

0 

0 " 


“o' 

0 

0 

0 

0 

1 
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0 

0 

0 

0 

0 

0 

1 
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0 
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0 

0 

0 

2n 

0 


Wx 

0 

0 

0 

-2n 

0 

0 


Uy 

0 

0 

-n^ 

0 

0 

0 J 


_“z_ 


(3.9) 


(3.10) 


Two of the primary assumptions of this model are a spherieal earth with no other 
perturbing forces present and a circular reference orbit. Both of these assumptions lead 
to unstable formations when applied to actual orbits since their effeets are signifieant. 
The latter assumption is the primary focus of this thesis and will be addressed in a later 
chapter. The third assumption is 

r«i? (3.11) 


whieh remains valid for most formations even when applied to aetual orbits. 

C. C-W SOLUTIONS 

For all solutions shown in this thesis, the filled cireles show the node points 
corresponding to the solution. They vary in number based on an arbitrary user-defined 

specification. The nature of the solver is sueh that the spacing between node points is not 
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constant. Instead, the points are closer together at the endpoints than they are in the 
middle. The solid lines represent the motion of the propagated initial eonditions. The 
solutions for the C-W equations of motion are not shown for the sake of brevity since 
they are identical, within numeric tolerances, to the solutions for the new equations 
shown later. 

1, C-W Circular Solution 

For the C-W equations, imposing the eondition of full periodicity and restricting 
the inter-satellite range to a constant distance from the eenter X (at point [0,0,0]) 
produees a cireular satellite formation. 

The final time for the optimal period was not fixed. The solution was not only 
able to determine the optimal trajectory, but also the optimal period in which to complete 
its trajectory. As expected, the solution resulted in = 1.0, meaning the optimal period 

was exactly equal to one orbital period. Interestingly, when asked to find a solution over 
2 orbits, the same 1-orbit solution was found, only it was now shown over the 2 orbits. 

2, C-W Projected Circular Solution 

Another well-known solution to the C-W equations is the projected circular 
formation. This solution maintains the appearance of a eircular formation as seen from 
the surface of the earth, but is in faet elliptieal. To aecomplish this formation a different 
path eonstraint, g, was imposed. 

(r-5)^ < -l-r/) < (r-l-5)^ (3.12) 

It is no longer the three dimensional range that was constrained but the range from 
the reference point to the satellite in a certain plane, namely the cross-track versus along- 
traek plane. As with the eircular formation, the optimal period was not fixed, and again it 
was exaetly equal to the orbital period. 

D. NON-LINEAR TWIST 

To further amplify the notion that linear models are not necessary for this 
approach. Thrust was chosen as a control variable instead of aeceleration. Not only is 
this more realistic but it also makes the “linear” equations “nonlinear” due to a non-zero 
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mass flow rate. This requires the addition of Mass as one of the state variables, in 
addition to the relative position and velocity components. 


X = 




(3.13) 


The path constraint remain s the same as equations (3.6) or (3.12). The state 
constraints now include the following constraint on mass: 

0.10<m(T)<1.0 (3.14) 


which assumes a minimum mass of 10 percent of the original mass and a maximum mass 
equal to the original mass. The event constraints also contain the new event 

m(To) = 1.0 (3.15) 

to force a starting mass at the beginning of the solution. 

The controls represent thrust in a given axis and are constrained by 

7;,r,,r.<r.„ (3.i6) 


with usually chosen in such a way as to make the constraint inactive. The desired 
cost is the total amount of thrust expended and can be described as 

J = - \%\ + \T\ + \T\dx (3.17) 

but was implemented as 

J = - \[T^ + Ty+T^)dx (3.18) 

This cost function could have been implemented as seen in equation (3.17) but, for 
numerical performance reasons, equation (3.18) was chosen. If the cost in equation 
(3.17) is minimized, the formation solution will be identical to the formation found by 
minimizing the cost in equation (3.18), when T.=0. For zero-thrust solutions, the main 
difference is the speed of the process. 
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1, New Equations of Motion 

The equations of motion for this ease are similar to Equation (3.8) except the 
control accelerations, u, are replaced with thrust. 



■ 0 

2n 
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"3n" 
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(3.19) 


Another equation must be included to govern the new state variable of mass. This mass 
flow rate is 


-T 



(3.20) 


Note that T is not the magnitude of the thrust vector, but is defined as 
T = |T)| + |T|,| + |t^ I and is the characteristic exhaust velocity of the thruster. 


2. Circular Formation Solution with New EOM 

The circular formation solutions to the new EOM are obtained using the path 
constraint shown in equation (3.6). Table III-l shows the constraints used for this 
solution. Eigure III-l shows the three-dimensional view, while Eigure III-2 to Eigure 
III-4 show the orthogonal projections. Eigure III-3 also shows the planar motion in the 
x-z plane at an angle of 60° to greater than 10'^ accuracy. Eigure III-5 shows the thrust 
profile, which is equal to zero-thrust within numerical tolerances. 


20 



Cross-Track 


Table III-l Constraints for Circular Formation 


Constraint 

Normalized 

Lower and Upper Bounds 

States: C’C’C 

Unconstrained 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

Unconstrained 

Time: x 

Unconstrained 

Events: r(Xo)-r(x^) 

[0.0 : 0.0] 

Events: f(Xo)-f(X|^) 

[0.0 : 0.0] 

Path: g = r^ + + r/ 

[1.0 : 1.0] 

Number of Nodes 

120 

Reference Orbit: e 

0.0 



Figure III- 1 Circular Formation Using New Equations of Motion 


21 










Cross-Track 
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Along-Track 


Figure III-2 Radial vs. Along-Track Motion for a Circular Formation 



Radial 


Figure III-3 Cross-Track vs. Radial Motion for a Circular Formation 




Figure III-5 Thrust Profile for Cireular Formation with New EOM 
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3, Projected Circular Formation Solution with New EOM 

The projected circular formation solutions utilize the new EOM with the path 
constraint in equation (3.12). Table III-2 shows the constraints while Figure III-6 shows 
the three-dimensional view. Figure III-7 to Figure III-9 show the orthogonal projections 
and Figure III-IO shows the thrust profile. 


Table III-2 Constraints for Projected Circular Formation 


Constraint 

Normalized 

Lower and Upper Bounds 

States: r^,ry,r^ 

Unconstrained 

States: r^,ry,r. 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

Unconstrained 

Time: x 

Unconstrained 

Events: r(xQ)-r(x^) 

[0.0 : 0.0] 

Events: f(xQ)-f(x^) 

[0.0 : 0.0] 

Path: g = ry + r/ 

[1.0 : 1.0] 

Number of Nodes 

100 

Reference Orbit: e 

0.0 
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Radial 


Figure III-6 Projected Circular Formation Using New Equations of Motion 



Figure III-7 Radial vs. Along-Track Motion for a Projected Circular Formation 












Figure III-IO Thrust Profile for Projeeted Cireular Formation with New EOM 


E, VALIDATION OF SOLUTIONS 

Sinee the answer to these problems was known, the validation was simply to 
compare the numerical solution to the analytical solution. This was accomplished using 
Figure III-l through Figure III-10 as well as unique characteristics such as the angle of 
the plane of relative motion seen in Figure III-3. This angle was calculated for each 
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solution and verified against the known values. 

Validation was also accomplished by numerically propagating the initial 
conditions for a set number of orbits, arbitrarily chosen as 50 orbits and is shown in 
Figure III-l 1 and Figure III-12. Table III-3 shows the results from propagating the initial 
conditions of the analytical solution. Table III-4 and Table III-5 show the results for 
propagating the initial conditions from the DIDO solution for the circular formation and 
projected circular formation, respectively. The errors seen in Table III-3 form the 
baseline for interpreting the errors in every other solution propagation. The source of 
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these errors is the numerie propagator, since the equations of motion used to propagate 
are the same as the analytical model. 

As mentioned previously, the optimal period was not fixed for these problems. 
Each of the solutions in this chapter yielded an optimal period equal to the orbital period, 
which was further verification of agreement with the known solutions. 
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Cross-Track 



Figure III-12 Projected Circular Solution Over 50 Orbits 


Table III-3 Propagation Results for Analytical C-W Circular Solution 


State Variable 

% Difference Between 
Initial and Final Values 

C 

9.84 xlO'^ 

c 

1.33 xlO'^ 

c 

9.84 xlO'^ 

K 

1.33 xlO'^ 

c 

9.84 xlO'^ 

c 

1.33 xlO'^ 
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Table III-4 Propagation Results for Cireular Formation 


State Variable 

% Difference Between 
Initial and Final Values 

C 

9.32 xlO'^ 

c 

3.15 xlO'^ 

c 

9.32 xlO'^ 

c 

3.15 xlO'^ 

c 

9.32 xlO'^ 

K 

3.15 xlO'^ 


Table III-5 Propagation Results for Projected Circular Formation 


State Variable 

% Difference Between 
Initial and Final Values 

C 

9.30x10'^ 

c 

3.13 xlO'^ 

c 

9.30x10'^ 

c 

3.16x10'^ 

c 

9.30x10'^ 

K 

3.16x10'^ 
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IV. ELLIPTIC REFERENCE ORBIT 


A. PROBLEM DEFINITION 

In many applications, it is desirable to design formations with non-cireular 
referenee orbits. The equations of motion from Chapter III are inapplieable for elliptical 
reference orbits.'^’ This chapter addresses these formations; in particular, the design of 
new formations using the method deseribed in Chapter IT 

Using the same eoordinate system as the previous ehapter, an additional state was 
added to simplify the equations of motion. Adding v , or True Anomaly, the state vector 
becomes 



The eontrols, representing both a positive and negative thrust direetion for eaeh 
axis, are 

" = [7’„ r,, T^_ 7;_]' (4.2) 


and were eonstrained by. 




T T T 

’ y+’U-’ ^ 


z+’ 


T. <T 


(4.3) 


The rationale for ehoosing the eontrols above was mainly for mathematical and 
numerical purposes. In ealeulating both the mass flow and eost of the solution, it is 
neeessary to ealeulate a total thrust value. Normally this is done by taking the sum of the 
absolute value for the individual thrust variable whieh presents mathematieal challenges 
due to non-differentiability. That is, the absolute value of any variable is defined as 


X 


A 


X 

-X 


x>0 

X < 0 


(4.4) 


It is evident that at x = 0 , the absolute value function becomes nondifferentiable. 
This ean be avoided by using 


X 
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(4.5) 



but even this introduees numerieal issues sinee the derivative of the square root function 
is undefined at x = 0, which is of great concern since = 0 is not only possible, but 
desired. 

One workable solution is to use a different thrust variable for the positive and 
negative directions. This choice models real life more closely since actual thrusters are 
capable of firing only in one direction. More importantly, it removes the numerical 
difficulties by allowing total thrust to be defined as 

r = I u'- = 7 ;,+ T^_+ 7 ;,+7;_+ 7 ;,+7;_ (4.6) 

This implementation does produce an additional concern since 71 = 0 is no longer 

internal to the constraints on the controls therefore making the constraints active. 
Activation of the Lagrange multiplier for the controls does not change the problem, but it 
does make the analysis and validation a bit more involved. Note that T does not 
represent the magnitude of the thrust vector. It has been redefined according to 
equation (4.6). 

For this formation, the path constraint is defined by, 

+ + Vx (4.7) 

where, di and du are minimum and maximum distance from the reference point to any 
spacecraft in the swarm. The cost function is 

1 

J = - \TdT (4.8) 

v-'o,. 

It is obvious that if the optimal cost turns out to be zero, then the solution corresponds to 
a zero-propellant formation; otherwise, the optimal (i.e. minimum fuel) open-loop 
controls to achieve the desired formation are obtained. 

B, EQUATIONS OF MOTION 

Since dynamics from Chapter III are invalid for an elliptical reference orbit, a new 
set of equations must be used. These equations of motion are described in References 17 


32 



and 27 and are derived here for the purpose of eompleteness. Following Kane’s 
notation, the superseript N is used to denote the Newtonian or Inertial reference frame. 
Starting with the most general equation describing the motion of the unperturbed 
reference orbit, 


K=-^ 


iRf 


(4.9) 


The motion of a swarm satellite is given by 


R:„=Rr.,+f''=-n7^+»-+»" 


R + r 


p T 


(4.10) 


where 


R + rf = (|R|'+2(R.r)|R| + |r|"|R|)' 


(4.11) 


If the assumption R » r is now introduced. 


R + r I 


\R + rf IRP 


„ ^ Rt „ 

R + r — 3 —^R 

IrI 


(4.12) 


The relative dynamics, though still in the inertial frame, can now be described by 




1 ^ 


R 


IRI' 


+ K+^t 


(4.13) 


J 


In order to retrieve the relative dynamics in the formation reference frame, the 
Transport Theorem for moving coordinate systems must be applied. The formation 
reference frame is denoted by the superscript B . Generally it is 


=f^ + (*'e)^xr) + (^e)^x'^e)^xr) + 2(*'e>^xf^) (4.14) 

Solving for f ^ yields 

f^=f^-(Wxr)-(Wx^e>^xr)-2(^e>^xf^) (4.15) 
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For an unperturbed orbit around a eentral body, using the perifocal inertial coordinate 
system, 


w 


[0 0 v] and '^e>''=[0 0 


• iT 


(4.16) 


The above equations imply that there cannot be a satellite at the reference point, since it 
is defined as an unperturbed reference orbit. The definition of angular momentum for 

90 

any unperturbed orbit. 


h = R^v = Vl 


— e 


solving for v , 


na 


v = —:-^/l-e^ with a 


i?(l + ecos(v)) 

( 1-4 


finally yielding 


n(l-ecos(v))^ 


(i-’f 

Differentiating equation (4.19) with respect to time produces 

.. -2nesin(v)v 


V • 


(- 4 ' 


-(l + ecos(v)) 


One additional term, , is defined as follows 


(4.17) 


(4.18) 


(4.19) 
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Substituting equations (4.13), and (4.16) through (4.21) into equation (4.15) 
provides the relative dynamieal equations of motion in the relative coordinate frame. 
They are 


with 



■ 0 

2v 

0" 


"v'+2^ V 

0 " 



f" = 

-2v 

0 

0 

f"+ 

-V v^-k 

0 

1 

r-l- — 
m 

X 

Ty 


0 

0 

OJ 


1 

o 

o 

-k 




(4.22) 


k = ri 


l + ecos(v) 

Tv 


(4.23) 


and an additional dynamic equation needed for the mass flow, 

m = — (4.24) 


One significant fact is that these equations are non-linear in mass. In order to 
limit the scope of this chapter, problems were confined to a zero perturbation case by 
setting = 0 . It should be noted, however, that the effects of any modelable perturbing 

force may be included in as desired to increase the fidelity of the model. 


C. NATURAL SOLUTIONS 

The first set of solutions are Natural solutions for varying eccentricities of the 
reference orbit. Natural is meant to denote that they require no thrust to maintain 
configuration and are at least partially periodic satisfying the traditional formation 
preconceptions. They are of course, only as good as the dynamical model used to find 
them. 


1. Natural Formation for e = 0.5 

The solution was found in normalized units, allowing it to be applied for any 
desired formation distance, RU. For the purpose of clarity and comparison, the 
solutions will be presented in nondimensional units. Table lV-1 itemizes the constraints 
in place for this formation. 
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Figure IV-1 shows the three-dimensional view of the formation relative motion. 
Figure IV-2 shows the orthogonal projection of the relative orbit in the radial versus 
along-track plane. Figure IV-3 shows the projection in the cross-track versus radial plane 
while Figure IV-4 projects the orbit onto the cross-track versus along-track plane. As 
with the previous figures, the X is the location of the reference point. From these results, 
position appears to be naturally constrained. In other words, the relative orbit is fully 
periodic even though the constraints are imposed for partial periodicity in velocity. This 
can be explained by the fact that a fully periodic solution would satisfy the constraints for 
a partially periodic solution. Additionally, since velocity was constrained, natural orbital 
motion will tend to drive position toward the periodic solution. This is not a law, since it 
is possible to find a formation that is periodic in velocity but not in position. On the other 
hand, if the position is constrained with velocity free, the formations will return to their 
original position but the velocity may be so large as to make this return impossible. This 
is of course, “undesirable” for formation configurations. 


Table IV-1 Constraints for Natural Formation with e = 0.5 


Constraint 

Normalized 

Lower and Upper Bounds 

States: C’C’C 

Unconstrained 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

[0 : Unconstrained] 

Time: x 

[0:5] 

Events: r(Xo)-r(x^) 

Unconstrained 

Events: f(xQ)-f(x^) 

[0.0 : 0.0] 

Path: g = r^+rl+r^ 

[2 : 22] 

Number of Nodes 

199 

Reference Orbit: e 

0.5 
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Figure IV-1 3-Dimensional Formation Trajectory for e = 0.5 Natural Formation 
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Figure IV-2 Radial vs. Along-Track Motion for e = 0.5 Natural Formation 
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Figure IV-5 Thrust Profile for e = 0.5 Natural Formation 

Of eourse, all these results are within numerical tolerances arbitrarily chosen to be 
10' in normalized units. It turns out that the error in the position vector is within 10' 
units after one orbit. From Figure IV-5, above, it is clear that this is a zero-propellant 
formation configuration. 

2. Natural Formation for e = 0.3 

Solution Two is another new formation configuration using the constraints shown 
in Table IV-2, specifically with e = 0.3. Similar to solution one, the time unit, TV , was 
calculated for a reference orbit with a perigee altitude equal to 1000 km. One of the 
differences between this formation and the previous, aside from reference orbit 
eccentricity, is the periodicity constraints. This solution was constrained to be periodic in 
position with velocity free. The previous solution, in section 1, was periodic in velocity 
with position free. Figure IV-6 shows the three-dimensional trajectory for this orbit in 
the formation reference frame. Figure IV-7 to Figure IV-9 show the orthogonal 
projections of the relative orbit. Again, the X in the figures is the location of the 
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reference point. Finally, Figure IV-10 demonstrates that this is another zero-propellant 
formation. 
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Table IV-2 Constraints for Natural Formation with e = 0.3 


Constraint 

Normalized 

Lower and Upper Bounds 

States: C’C’C 

Unconstrained 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

[0 : Unconstrained] 

Time: x 

[0:5] 

Events: r(Xo)-r(x^) 

[0.0 : 0.0] 

Events: f(Xo)-f(x^) 

Unconstrained 

Path: g 

Unconstrained 

Number of Nodes 

80 

Reference Orbit: e 

0.3 



Figure IV-6 3-Dimensional Formation Trajectory for e = 0.3 Natural Formation 
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3. Natural Formation for e = 0.7 

Recently, Inhalan et al presented new periodic formations by an analytic method. 
Solution Three, shown in Figure IV-11 and Figure IV-12, reproduces a particular 
formation given in Reference 17 for e = 0.7 using this method. 

By finding the same solution, the method is again validated against independent 
results. The solution was found by constraining the position states to a three-dimensional 
box slightly larger than the proposed solution and using the fully periodic event 
constraints (see Table IV-3). 


Table IV-3 Constraints for Natural Formation with e = 0.7 


Constraint 

Normalized 

Lower and Upper Bounds 

States: 

[±0.62, ±1.4, -5.8:1.0] 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

[0 : Unconstrained] 

Time: x 

[0 : 10] 

Events: r(xQ)-r(x^) 

[0.0 : 0.0] 

Events: f(Xo)-f(x^) 

[0.0 : 0.0] 

Path: g 

Unconstrained 

Number of Nodes 

99 

Reference Orbit: e 

0.7 
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D, FORCED SOLUTIONS 

This section provides a glimpse into the power of this method. Specifically, that 
it can find not only natural formations, but also controlled or forced formations. These 
solutions are unique in the sense that Keplerian orbital dynamics cannot produce these 
formations since they must be actively controlled. This method will determine the 
feasibility of the formation and will automatically provide the open loop controls required 
to maintain the desired configuration. 

1, Forced Circular Formation 

The formation in this section uses a reference orbit eccentricity of 0.3. The 
configuration constraint is a fixed radius from the reference point. That is, the swarm is 
restricted to a surface that lies on the sphere r = constant. This path constraint is written 

g = (fi(4-25) 

The goal is to minimize the fuel required to meet the configuration constraints 
shown in Table IV-4. Figure IV-13 shows the three-dimensional plot of the solution 
formation. It closely resembles Figure III-l, the C-W circular formation, but closer 
inspection will reveal the subtle differences. Figure IV-14 through Figure IV-16 show 
the projection of the formation in the three orthogonal planes. Figure IV-17 shows the 
open loop controls required by one of the satellites to maintain this formation for the 
given eccentricity. This solution was found using 100 nodes and fixing the final time to 
exactly one orbit. 
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Table IV-4 Constraints for Forced Circular Formation 


Constraint 

Normalized 

Lower and Upper Bounds 

States: C’C’C 

Unconstrained 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

[0 : Unconstrained] 

Time: x 

x^ =1.0 

Events: r(Xo)-r(x^) 

[0.0 : 0.0] 

Events: f(Xo)-f(x^) 

[0.0 : 0.0] 

Path: g = r^ + + r/ 

[1.0 : 1.0] 

Number of Nodes 

100 

Reference Orbit: e 

0.3 
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Cross-Track ... Radial 



ure IV-14 Forced Circular Formation Radial vs. Along Track Motion 
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Figure IV-15 Forced Circular Formation Cross-Track vs. Radial Motion 
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A comparison of Figure IV-15 and Figure III-3 shows the departure from the 
classic C-W solution. Figure IV-15 apparently shows non-planar motion, but itself is not 
enough to determine if the formation actually resides in a plane. The answer lies in the 
angular momentum veetor, =rxf . The direction of the relative angular moment 

veetor was plotted to determine if the motion was planar. If indeed the motion is planar, 
then the direetion of the angular momentum should remain constant. If the formation is 
not planar, then the direction of the momentum vector will not be constant. Figure IV-18 
shows this plot for the classic C-W circular formation described in section III.D.2 and 
demonstrates that it is constant, within numeric tolerances. On the other hand. Figure 
IV-19 shows the same traee for this formation. It is clearly not constant and therefore 
demonstrates that this formation is non-planar. In both figures, a line is drawn from the 
origin in the direction of the unit vector associated with the relative angular momentum at 
the first time step. The dots depict the tip of this angular momentum unit veetor at each 
time step. 



Figure IV-18 Relative Angular Momentum Veetor for the C-W Circular Formation 
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Figure IV-19 Relative Angular Momentum Vector for the Forced Circular Formation 


2, Forced Projected Circular Formation 

This formation in another example of a forced formation using a well-known 
configuration. The configuration constraint is the same as seen in section III.D.3, except 
for a reference orbit eccentricity of 0.3. The goal remains the same: to minimize the fuel 
required to meet the configuration constraints shown in Table IV-5. Figure IV-20 shows 
the three-dimensional plot of the solution formation. Figure IV-21 through Figure IV-23 
show the projection of this formation in the three orthogonal planes. Figure IV-24 shows 
the open loop controls required by one of the satellites to maintain this formation for the 
given eccentricity. This solution was also found using 100 nodes and fixing the final 
time to exactly one orbit. 
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Table IV-5 Constraints for Forced Projected Circular Formation 


Constraint 

Normalized 

Lower and Upper Bounds 

States: C’C’C 

Unconstrained 

States: 

Unconstrained 

States: m 

[0.1 : 1.0] 

Controls: T 

[0 : Unconstrained] 

Time: x 

x^ =1.0 

Events: r(xQ)-r(x^) 

[0.0 : 0.0] 

Events: f(Xo)-f(x^) 

[0.0 : 0.0] 

Path: g = ry + r/ 

[1.0 : 1.0] 

Number of Nodes 

100 

Reference Orbit: e 

0.3 
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Figure IV-20 Forced Projected Circular Formation for e = 0.3 
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Figure IV-21 Forced Projected Circular Formation Radial vs. Along Track Motion 
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Figure IV-22 Forced Projected Circular Formation Cross-Track vs. Radial Motion 








E. VALIDATION OF SOLUTIONS 

For all the cases presented here, the following steps were taken in validating the 
solutions. The Lagrangian was constructed in accordance with section II.B.l and 
equation (2.19). Care must be taken to include the state and control bounds in defining 
the path constraints, g . Taking the partial derivatives of the Lagrangian with respect to 
the controls yields the following equations: 


ac 

1 

A 


(4.26) 


T/ — Tq 

m 

dL 

1 


- — + ^T 

(4.27) 

K- 

T j- — Tq 

m 

dL 

1 

K 
-1_ 

Ve 

(4.28) 


Xf — Xq 

m 

dL 

1 



(4.29) 

K- 

Xj- — Xq 

m 

dL 

1 



(4.30) 


Xf — Xq 

n- 

m 

dL 

1 


- — + ^T 

(4.31) 

dT,- 

Xf — Xq 

m 


(|) is called the switching function for the controls and is governed by the KKT conditions 
described in equation (2.22). Figure IV-25 shows the switching function for the Thrust in 
each of the three axes corresponding to the thrust profile shown in Figure IV-24. Not 
inherently obvious is that the switching function is the same for both circular and 
projected circular formations at given values of v, since the path constraints are not 
control dependent. 
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Figure IV-25 Switching Functions for Control Thrust 


In addition to the optimality analysis described above, the initial conditions from 
the DIDO solution were numerically propagated for 50 orbits. For the natural solutions, 
the initial conditions were propagated with all components of thrust set equal to zero. For 
the forced solutions, the solution thrust profile was imported into the propagator. This 
profile was interpolated to determine the value for the controls at each propagation time 
step. The propagator was free to determine its own time steps, which did not necessarily 
match the solution time steps. In the process of creating a propagator, every MATLAB 
resident ODE solver was evaluated. ODE45 proved to have the best combination of 
accuracy and speed. Additionally, many different interpolation schemes were evaluated. 
A MATEAB function called POEINT provided the most accurate results. It was created 
by Weideman and Reddy and implements the barycentric formula from Henrici’s 
Essentials of Numeric Analysis. Spline interpolation was the next best in accuracy, and 
provided a significant benefit in speed over POEINT. Spline was used for routine 
evaluations and POEINT was reserved for more detailed assessments. 
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1. Natural Formation for e = 0.5 

Figure IV-26 and Table IV-6 show that there is no appreeiable deviation in 
formation configuration after being propagated for 50 orbits. 



Figure IV-26 Natural Formation with e = 0.5 Over 50 Orbits 


Table IV-6 Propagation Results for Natural Formation with e = 0.5 


State Variable 

% Difference Between 
Initial and Final Values 


8.54x10'^ 

r 

y 

5.77 xlO'* 


4.61 xlO'^ 


6.67 xlO'* 


1.04x10'^ 

K 

2.35 xlO'^ 
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2, Natural Formation for e = 0.3 

Again, the initial conditions were numerically propagated (with thrust equal to 
zero) for 50 orbits. Figure IV-27 and Table IV-7 show that there is no significant 
difference over this period. 



Figure IV-27 Natural Formation with e = 0.3 Over 50 Orbits 


Table IV-7 Propagation Results for Natural Formation with e = 0.3 


State Variable 

% Difference Between 
Initial and Final Values 


9.06 xlO'^ 

r 

V 

8.71 xlO'^ 


1.97x10'^ 

K 

5.73 xlO'^ 


1.11 xlO'^ 

K 

6.67 xlO'^ 
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3. Natural Formation for e = 0.7 

Since solution three is a reproduction of a previously known solution, the 
validation of the formation configuration can be found in Reference 17. 

4, Forced Circular Formation 

The validation of this result included propagating the initial eonditions 
numerically for 50 orbits, now subject to the thrust profile shown in Figure IV-17. Figure 
IV-28 shows the three-dimensional motion of the formation over this time period. Table 
IV-8 also demonstrates that, as desired, there is no appreciable digression in the 
formation configuration over 50 orbits. 
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Figure IV-28 Foreed Cireular Formation over 50 Orbits 
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Table IV-8 Propagation Results for Forced Circular Formation 


State Variable 

% Difference Between 
Initial and Final Values 


0.23 


1.51 

C 

0.04 

K 

0.22 


0.80 

K 

0.11 


One unique consequence of imposing the circular formation constraint is that it 
provides another way to verify the results. The forced circular is defined by r = constant 

which means f = 0 . This implies = r ■ r = const, which gives ) = 2(f-r) = 0, 

followed naturally by 

— (f-r) = f-r + ff = 0 (4.32) 

dt 

Substituting the EOM for f above provides an independent check that the solution is 
indeed a circular formation. 

5. Forced Projected Circular Formation 

As with the previous solutions, the initial conditions were propagated numerically 
for 50 orbits, subject to the thrust profile shown in Figure IV-24. Figure IV-29 shows the 
three-dimensional motion of the formation over this time. Table IV-9 also details the 
errors in the formation configuration after 50 orbits. 
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Figure IV-29 Forced Projected Circular Formation over 50 orbits 


Table IV-9 Propagation Results for Forced Projected Circular Formation 


State Variable 

% Difference Between 
Initial and Final Values 

C 

0.06 

c 

1.33 

c 

0.01 

c 

0.11 

c 

0.47 

c 

0.02 
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V. UNRESOLVED ISSUES 


A. OPTIMAL PERIOD VS. ORBITAL PERIOD 

For Natural or Zero-Thrust solutions, the optimal period, is exaetly 

equal to the orbital period. However, onee thrust is required to maintain the formation 
the optimal period may or may not equal the orbital period. One way to resolve this is to 
foree the solution period to be equal to the orbital period. However, this fix raises several 
questions. 

Table V-1 shows the eost assoeiated with two different solutions to the Foreed 
Cireular Formation with e = 0.3 and assumes a 100 kg satellite, with 1^^ = 1000 see , and 

the TV for a perigee altitude of 1000 km. Solution lA is detailed in seetion IV.D.l, 
while Solution 2A is not shown in the prior ehapters. Table V-2 shows the eost 
assoeiated with two different solutions to the Foreed Projeeted Cireular Formation using 
the same assumptions deseribed above. Solution IB is represented in seetion IV.D.2, 
while Solution 2B is not shown. 


Table V-1 Differing Costs for Foreed Cireular Solutions 


Forced Circular Solution lA 

Cost: 

1.98 xlO'^ Newtons/see 

Cost: % Mass 

2.34x10'^% 

Optimal Period: 

1.0 orbital periods 

Forced Circular Solution 2A 

o 

o 

VI 

1 

1.80 xlO'^ Newtons/see 

Cost: % Mass 

1.00x10'^% 

Optimal Period: x^ 

1.0005 orbital periods 
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Table V-2 Differing Costs for Foreed Projeeted Cireuiar Solutions 


Forced Projected Circular Solution IB 

Cost: 

2.87 xlO'^ Newtons/sec 

Cost: % Mass 

3.39x10'^% 

Optimal Period: 

1.0 orbital periods 

Forced Projected Circular Solution 2B 

Cost: 

2.85 xlO'^ Newtons/sec 

Cost: % Mass 

3.38 xlO'^ % 

Optimal Period: x^ 

1.002 orbital periods 


In both oases, the seoond solution has the lower total oost. However, the optimal 
period of these solutions is not equal to the orbital period. The preoise explanation is 
unknown. These solutions may be periodio as a formation, but not periodio with respeot 
to the reference point since the period of the reference point (or satellite) is, by definition, 
equal to the orbital period. This mismatch in periods presents difficulty in visualizing or 
plotting the configuration or relative motion. If ignored, the formation appears to drift 
away from the reference point and therefore not stay together. One way to overcome this 
visualization problem is to display the position of one swarm satellite against another 
swarm satellite. Over the course of the optimal period, the distance should remain 
bounded and repeat itself 

B, SWITCHING FUNCTION : DERIVED VS. DIDO SOLUTION 

Normally, (see equations (4.26) through (4.31) ) is equal to zero for an 

dT 

optimal solution. This allows the values for the switching function, (|), to be derived and 
calculated. At the same time, DIDO calculates the values for the switching function 
directly. Figure IV-25 shows the switching function as calculated by DIDO, which 

agrees with the expected results. The concern is that the value for is not 0.0 but 

dl 
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exactly 0.5. The resulting switching functions have the identical shape as those in Figure 
IV-25, but are offset by negative 0.5 in the vertical axis. This discrepancy in the 
switching function requires further analysis. 

C. J2 PERTURBATIONS 

Work was started on implementing the effects of earth oblateness or J 2 in the 
relative frame, but due to time constraints was not completed. 

1, Linear J 2 terms. 

The following Linear equations for J 2 effects in the relative frame came from an 
unpublished paper by I.M. Ross. Using the same coordinate system as Figure II-l and 
assuming a circular reference orbit, the J2 effects can be represented as relative 
accelerations. 



12sm^ zsin^ 05 

-4 sin^ i sin 20 

-4sm2zsm0 

^7 

^^2 ^ R 

-4 sin^ i sin 205 

1-1-sin^ z (2 - 7 sin^ 0 ) 

sm2zcos0 



-4 sin 2zsm0 

sin 2z cos 0 

3-sm^ z(2-l-5sm^ 0 ) 



(5.1) 


where 


J. 


3/2^; 


and n5 = co+v 


(5.2) 


Non-linear J 2 terms. 


20 

By comparison, the non-linear J 2 terms in the inertial reference frame are 
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where 


J. 


2 


(5.6) 


In order to get the relative aecelerations resulting from J 2 into the relative frame, several 
transformations must be defined. To transform perifocal to Earth-centered inertial (ECI): 


T Tjy 


cQcco- sQs(Oci 

-cQscd- sQccoci 

sQsi 

UK 

— 

sQco)+cQs(Oci 

sQs(£)+cQc(£)ci 

-cQsi 

PQW ' 







scdsi 

ccdsi 

ci 


with c and s representing sin and cosine of the appropriate angles. The transpose of 
above is used to convert from ECI to perifocal. Transforming the formation reference 
frame to perifocal is done by 




cosv 

-sinv 

0“ 

~ PQW\ 




0 

_RSW \ 

— 

sinv 

cosv 



0 

0 

1 


with the transpose used for moving in the other direction. 


The algorithm used to calculate J 2 from the relative positions of the swarm 
satellites begins with converting relative positions into inertial positions. 


~ PQW\ 


~PQW~ 

~ RSW~ 

UK ) 

UK 

PQW 


(5.9) 


The next step is to calculate the accelerations in the inertial frame according to equations 
(5.3) through (5.5). These inertial accelerations can now be transformed, using equation 
(4.15) into the relative accelerations. Of note is the use of a static reference orbit. That 
is, the reference orbit is assumed to be unperturbed. 

3, Comparison Between Sets of Equations. 

Initial results from a comparison between the linear equations and the non-linear 

equations were encouraging. The comparison was done for a circular formation, with 

e = 0 at three different inclinations: 28.5°, 45°, and 63.4°. Eigure V-1 to Eigure V-3 show 

the perturbation accelerations due to J 2 effects for all three inclinations, normalized by the 
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scaling scheme detailed in seetion II.B.2. Figure V-4 to Figure V-6 show the difference 
between equation (5.1) and equations (5.3) through (5.5) for each inclination. Note that 
the differences are not random, but seem to follow some pattern. This is the effect of the 
higher order terms not present in the linear equations. 



Figure V-1 Relative Aceelerations Due to J2 at 28.5° Inclination 
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Relative J„ Accelerations Relative Accelerations 
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Figure V-2 Relative Aeeelerations Due to J2 at 45“ Inelination 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 


Time 

Figure V-3 Relative Aeeelerations Due to J2 at 63.4“ Inelination 
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VI. CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

This work for this thesis started as a question: Can optimal periodic control 
(OPC) theory be applied to satellite formation design and control? Using the 
fundamentals of optimal control theory, a framework was developed that captures the 
essence of designing and controlling spacecraft formations. This framework is used to 
articulate a variety of formations including the notion of an aperiodic formation. Based 
on a deliberate formulation, including mass as a state variable, it was shown that the 
numerical approach can easily handle nonlinearities. Additionally, formations were 
presented that require active control along with their corresponding thrust profile. 

This thesis lends credence to the notion of numerically searching for minimum- 
fuel formation configurations for spacecraft swarm subject to arbitrary nonlinear 
dynamics. Thus, practical formations may be designed and controlled using this method. 
The foundation for applying OPC to satellite formations has been completed. Future 
work can build on this foundation in increasing the complexity of the models. A 
substantial amount of work needs to be done in developing a more complete model of 
satellite motion and perturbations. The framework developed here can be readily used 
with more robust dynamical models and should produce very interesting results. 

B, OPPORTUNITIES FOR FUTURE WORK 

In addition to the topics identified in the previous chapters, opportunities for 
future work are available in the following areas. 

In the course of implementation for Chapter IV, an additional “state” was added 
to simplify the equations. This state variable was v with its corresponding derivative, v . 
This required adding dynamical constraints for this state. One opportunity for future 
work would be to include v in the path constraints and rewrite the dynamical equations 
to remove explicit dependence on v. This may possibly increase numerical accuracy 
when implementing the model in DIDO. 


71 



Another possibility for follow on research is an examination of the ^6)^ term . It 
may be possible to continue working in the formation reference frame, while adding the 

J2 effects to this term by including — , and . 

dt dt dt 

A departure from the relative reference frame would be an area for future work. 
The creation of a model in an inertial frame and the use of the full non-linear equations 
will create a high fidelity model, and will allow the addition of J2 effects and any other 
Earth-centered effects to be included more readily. Translating the formation constraints 
from the relative frame, where they are intuitive, to the inertial frame should be 
straightforward but tedious. If the model is defined in inertial space, the deployment and 
reconfiguration problems becomes a simple change from one orbit to another orbit. The 
difficulty lies in specifying the constraints, if any, on the reconfiguration relative motion 
and visualization. 


Another consideration, offered by Dr. Terry Alfriend, is to write the path 
constraint as follows, and minimize total thrust. 




I J 


A 

V y 


-I- 


\R. 


'fef\ 


= constant 


( 6 . 1 ) 


If the reference orbit is circular, then is constant and becomes a simple scaling 

factor, which was completed previously. The above constraint was not addressed for an 
elliptical reference orbit, where is no longer a constant. 

Other unexplored areas include the possibility of near-circular relative formations. 
If a small deviation is allowed in the path, it may lower the cost of the formation. In 
addition, specific missions require specific configuration constraints. One very 
interesting area of research would be to identify and catalog various configuration 
constraints according to mission. 
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